Non-equilibrium Gross-Pitaevskii dynamics of boson lattice models 



(N 
O 
O 



43 

o 



CO 



I 

o 
o 



(N 
> 
O 
On 
^t- 
VO 
O 
(N 
O 



o 



X 
S3 



Anatoli Polkovnikov,[] Subir SachdevJ] and S. M. Girvin| 

Department of Physics, Yale University, P. 0. Box 208120, New Haven CT 06520-8120 

(Dated: February 1, 2008) 

Motivated by recent experiments on trapped ultra-cold bosonic atoms in an optical lattice poten- 
tial, we consider the non-equilibrium dynamic properties of such bosonic systems for a number of 
experimentally relevant situations. When the number of bosons per lattice site is large, there is a 
wide parameter regime where the effective boson interactions are strong, but the ground state re- 
mains a superfluid (and not a Mott insulator): we describe the conditions under which the dynamics 
in this regime can be described by a discrete Gross-Pitaevskii equation. We describe the evolution 
of the phase coherence after the system is initially prepared in a Mott insulating state, and then 
allowed to evolve after a sudden change in parameters places it in a regime with a superfluid ground 
state. We also consider initial conditions with a u n phase" imprint on a superfluid ground state 
(i.e. the initial phases of neighboring wells differ by n), and discuss the subsequent appearance of 
density wave order and "Schrodinger cat" states. 



I. INTRODUCTION 

With the emerging experimental studies of ultra-cold 
atoms inrSh parabolic trap and a periodic optical lattice 
potential!™ (the wavelength of the optical potential is 
much smaller than the dimensions of the trap), new possi- 
bilities for studying the physics of interacting bosons have 
emerged. At equilibrium, the bosons can undergo a tran- 
sition from a superfluid to an insulaljaL^Sjihe strength 
of the optical potential is increaseduu^EfLla. However, 
the facile tunability and long characteristic time scales 
of these systems also offer an opportunity to investigate 
non-equilibrium dynamical regimes that have not been 
accessible before. In this context, there have been a few 
recent theoretical studies of the dynamics of bosons in 
a periodic potential: Ref. ^| computed the oscillation 
frequency of the center of mass of a superfluid state 
of bosons, while, some-. non-equilibrium issues were ad- 
dressed in papersE30'E3 which appeared while this paper 
was being completed. 

A description of the purpose of this paper requires an 
understanding of the different parameter regimes of the 
boson system, which we will assume is well described by 
the single-band Hubbard model: 
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Here cij is a canonical Bose annihilation operator on sites 
of the optical lattice ("wells") labeled by the integer j, 
J is the tunneling amplitude between neighboring lattice 
sites, U > in the repulsive interaction energy between 
bosons in the same lattice minimum, and Vj is a smooth 
external potential which we will take to be parabolic. We 
will mainly consider the case of a one-dimensional optical 
lattice, relevant to the experiments of Ref. [j], but gener- 
alization to higher dimensions is possible. The form of 
Vj and the chemical potential of the bosons determine 



another important parameter: N, the mean number of 
bosons at the central site (more precisely, at the site 
where Vj is smallest); we shall mainly consider the case 
N 1 here. A dimensionless measure of the strength of 
the interactions between the bosons is the coupling 
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the different physical regimes of TL are also conveniently 
dilineated by the values of A. 

When the interactions between the bosons are strong 
enough, A > Asj, the ground state of TL undergoes a 
quantum phase transition from a superfluioLto a Mott 
insulator (see Appendix |A|). It is known thatl± 
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So for the case where N is large, there is a wide regime, 
1 <C A <C ./V 2 , where the interactions between the bosons 
are very strong, but the ground state is nevertheless a 
superfluid. A description of the dynamical properties of 
TL in this regime is one of central purposes of this paper. 

For iV— large, and A smaller than As/, it is widely 
acceptedty that the low temperature dynamics of TL can 
be described by treating the operator classical 
c-number. (We will investigate the conditions for the va- 
lidity of this classical approximation more carefully in 
Section O, where we will also discuss the time range 
over which it can be applied.) More precisely, we intro- 
duce the dimensionless complex dynamical variable ipj (t) 
whose value is a measure of (dj (t)) / yiV ; then its dynam- 
ics is described by the classical Hamiltonian 
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and the Poisson brackets 
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Here, and henceforth, we measure time in units of h/J. 
The resulting equations of motion are, of course, a dis- 
crete version of the familiar Gross-Pitaevskii (GP) equa- 
tions. We will often impose a parabolic confining poten- 
tial, in which case 
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A nonuniform potential Vj also can lead to localization 
of bosons in separate wells; in particular, even without 



where the 4>j are independent random variables which are 
uniformly distributed between and 2k. In this manner, 
we have mapped the fully deterministic quantum time 
evolution of Ti to the stochastic and classical time evolu- 
tion of Hqp. In practice, the procedure is then as follows: 
choose a large ensemble of initial values of 4>j > an d deter- 
ministically evolve Hqp for each such initial condition; 
the expectation value of any quantum observable at time 
t is then given by the average value of the correspond- 
ing classical observable at time t, with the average being 



interaction (A = 0), when \Vj + i 
modes of (1.4) become localized. Note that this local- 
ization is a purely semiclassical effect, described by the 
GP equations. If Vj is smooth then for A > As/, the 
system, undergoes a transition to nonuniform insulating 
stateOM 

Describing the non-equilibrium quantum Bosc dynam- 
ics for A < Xsi is now reduced to a problem of integrat- 
ing the classical equations of motion implied by ( |l.4|l-- r )| ) ■ 
However, it remains to specify the initial conditions for 
the classical equations; these clearly depend upon the 
physical situations of interest, and we shall consider here 
two distinct cases, which are discussed in the following 
subsections 



A. Mott insulating initial state 

ConsiderjJthe physical situation (of current experimen- 
tal interestEj) where for t < the bosons are in a Mott 
insulating state with A > Xsi , and at time t = the opti- 
cal lattice potential is suddenly reduced so that A < Xsi 
for all t > 0. Clearly, the GP equations should apply for 
t > 0, and the Mott insulating initial state will impose 
initial conditions which we now describe. The required 
initial conditions are readily deduced by thinking about 
the full quantum Heisenberg equations of motion for a,j (t) 
implied by Ti.. By integrating these equations, one can, in 
principle, relate any observable to the expectation values 
of products of powers of Oj(t = 0) and a,j(t = 0). For the 
Mott insulator with A 3> Xsi these expectation values 
have a very simple structure: they factorize into prod- 
ucts of expectation values on each site, and are non-zero 
only if the number of creation and annihilation operators 
on each site are equal. Furthermore, for large N, we can 
also ignore the ordering of the cij and a - operators on 
each site, and e.g. we obtain to leading order in 1/JV: 

(af{t = Q)a?(t = 0)} « S nrn 6 je (N 3 ) n , (1.6) 

where we have accounted for a possible spatial inhomo- 
geneity by introducing Nj (a number of order N), the 
number of bosons at site j in the Mott insulator. In 
terms of the classical variables tpj > the t — expectation 
values in (11.61) are easy to reproduce. We simply choose 



V \ > 2J the eigen- taken over the random variables <j>j. In particular 
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where we have indicated that the angular brackets on the 
left represent a traditional quantum expectation value, 
while those on the right represent a n av erage over the 
independent variables <pj specified by (1.7) at time t = 0. 
We will henceforth implicitly assume th at a ll angular 
brackets have the meaning specified in ([0]), depend- 
ing upon whether they contain quantum or classical vari- 
ables. 

An important property of (1.8) is that while we must 
have j ' = j for a non-zero result at t — 0, this is no longer 
true for t > 0. In particular, non-zero correlations can 
develop for large \ j ' —j \ as time evolves, corresponding to 
a restoration of phase coherence. Indeed the ground state 
for A < As/ is superfluid and thermalization must lead to 
increase of the phase correlations. However, in this paper 
we show, that even without relaxation the coherence can 
be restored dynamically. (Of course, as we are looking at 
one dimensional systems and the final state is expected 
to be thcrmalized at a non-zero temperature, the phase 
correlations cannot be truly long-range and must decay 
exponentially at large enough scales: however, guided 
by the experimental situation, we will look at relatively 
small systems for which this is not an issue.) Describing 
the dynamics of the restoration of this phase coherence is 
also a central purpose of this paper. We shall characterize 
the phase coherence by studying the expectation value of 
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where M is the number of lattice sites (for a nonuni- 
form external potential Vj, M is just the ratio of the 
total number of bosons to the number of bosons in cen- 
tral well), and g is some suitably chosen weight function. 
Observables closely related to T> g are measured upon de- 
tecting the atoms after releasing the trap. At time t = 0, 
T> g (0) = 0, and we will be interested in the deviations of 
T> g (t) from this value for t > 0, an increase correspond- 
ing to an enhancement of superfluid phase coherence. We 
note, iiLpassing, that a closely related procedure was used 
earlieiLj to describe the onset of phase coherence after a 
sudden quench from high temperature; here, we are al- 
ways at zero temperature, and move into a superfluid 
parameter regime by a sudden change in the value of A. 
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We will begin our analysis of the structure of T> g (i) 
by c onsid ering the case with two wells (M = 2) in Sec- 
tion II A. For the weakly interacting case (A <C 1), D g (t) 
exhibits Josephson oscillations with a period of order 
unity; the weak interactions lead to a decay of oscilla- 
tions with a slow (i -1 / 2 ) saturation of the coherence at 
a steady-state value at a time scale t oc A -1 . For A ^> 1 
the oscillations are suppressed and T) g (t) saturates at 
t oc 1/vA, which is, in fact, shorter than a single tun- 
neling time. For this two lattice site case we can also 
obtain a complete solution for T> a {t ) for the quantum 
Hamiltonian H (described in Section [I A 2), and this al- 
lows a detailed analysis on the regime of validity of the 
semiclassical GP equations. We show that the semiclassi- 
cal approach is valid for two lattice sites when N is large 
and t < N/X. This is, in fact, a general result which 
implies that the quantum mechanics becomes important 
when time exceeds inverse energy level spacing. For more 
than two lattice sites, the energy splitting scales as the 
inverse of the total number of particles and at A <C 1, the 
semiclassical conditions arc virtually always fulfilled. It 
is surprising that even with a small number of particles 
N = 4, and weak interactions, the GP equations give an 
excellent description of the system evolution, apart from 
overall numerical prefactor (1 + 2/JV), which is not small 
in this case. 

The restoration of coher ence is also studied in the 
many well case in Section III A. We discuss the case 



without an external potential in Section III A 1 ; with an 
equal number of particles initially in all the wells, phase 
correlations develop only in the interacting case (A > 0). 
This is true for both periodic and open boundary con- 
ditions. Similar to the two well case, in the weakly in- 
teracting regime phase correlations will oscillate in time. 
However these oscillations will be periodic only for partic- 
ular number of wells: M = 2, 3, 4, 6 for periodic boundary 
conditions and M — 2, 3, 5 for open boundary conditions. 
For other numbers of wells, the oscillations are chaotic. 
As for the two well case, a stronger interaction results in 
decay of correlations in time, leading to the steady state. 



Next, in Section III A 2, we consider the restoration of 



phase coherence for the experimentally important case of 
a parabolic potential. The results are quite different for 
this case, and phase correlations develop even without 
interactions. In a weak parabolic potential, T> g (t) oscil- 
lates with a frequency which scales as the square root of 
the parabolicity, £. This frequency is closely related to 
the. oscillation frequency discussed recently by Kramer et 
al.u for the case where the center of mass of the atomic 
gas is displaced. In the present situation, there is no 
displacement of the center of mass, but the same oscilla- 
tion is excited upon a sudden change in the value of A. 
The oscillations decay even at A = 0; weak or intermedi- 
ate interactions A < 1 do not change the noninteracting 
picture much. The amplitude of the oscillations become 
more pronounced for A « 1, but for A ^> 1 the oscillations 
are suppressed as for the flat potential. 

While this work was being completed, we became 



aware of related results of Altman and Auerbach also ad- 
dressing the restoration of phase coherence in a Mott in- 
sulator. However, there are some significant differences in 
the physical situations being addressed. Above, we have 
considered a system deep in the Mott insulating phase 
(with A ^> Agi) taken suddenly to parameters for which 
the ground state was deep in the superfluid phase (with 
A <C As/). In contrast, Ref. |l^ consider the case when 
both the initial and final values of A were not too far from 
As/, but remained on opposite sides of it. For A close to 
As/, and at temperatures not too small, a "relativistic 
Gross-Pitaevski" equation had been proposed in Ref. [l8| 
as a description of the "_Bo.se molasses" dynamics of the 
order parameter. The conditions under which oscillations 
in the amplitude of the ordet-parameter would be unde&. 
damped were also presented^. Altman and AuerbachliZl 
advocated that the same equations could describe the 
time evolution of the amplitude of the order parameter 
as it evolved from the Mott insulator (with zero ampli- 
tude) to the superfluid (with finite amplitude) at zero 
temperature. We review issues related to the damping of 
the amplitude mode in Appendix |B| Altman and Auer- 
baclJiJ also considered the situation without an external 
potential (Vj = 0). We have noted above that such a po- 
tential changed our results significantly; in Appendix ^ 
we discuss the significant role of the external potential in 
the equilibrium properties for A As/- 



B. Modulated phase initial state 

A second set of initial conditions we consider is the 
case in which the parameter values always correspond 
to a superfluid ground state i.e. A < As/- For time 
t < we imagine that A takes some fixed value and the 
phases <f>j have some known set of fixed, non-random val- 
ues at t = and we follow the subsequent evolution of 
the bosons using the discrete GP equation. The phase 
imprint can be experimentally achieved by e.g. applying 
a short (compared to a single tunneling time) pulse of 
external field to the condensate. A case of special inter- 
est will be when there is a relative 7r phase shift between 
neighboring wells: 
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(1.10) 



For two wells with equal Nj and relatively small A, this 
state is metastable (this is also the case for even M and 
periodic boundary conditions). However, if the interac- 
tion A becomes larger than a critical value, this equi- 
librium becomes unstabla.and the bosons spontaneously 
form a "dipole" stateE3E3E!J in which most of them oc- 
cupy one of the two wells (see Section II B). Upon ac- 



counting for quantum tunneling in a system with a finite 
number of bosons, the state obtained is a superposition 
of the two dipole states restoring translational symme- 
try. Howe ver, in case of infinite number of wells (see Sec- 
tion III Bj ) the tunneling between the two dipole configu- 
rations is negligible and translational symmetry is broken 
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by the appearance of a density wave of bosons with a pe- 
riod of two lattice spacings. This effect is similar to that 
studied in Ref. |2l] for the case of a Mott insulator in a 
strong electric field. 

Related to this instability is a very interesting pos- 
sibility of formi ng a Schrodinger cat stateEj. We show 
in Section [II B that if the system is initially in the "7r 
state" , and the interaction is slowly increased, then at 
certain point all the bosons spontaneously move into one 
of the wells. If quantum mechanical corrections are taken 
into account then the final configuration is the superpo- 
sition of the states with all bosons in one of the wells. 
This effect opens the possibility of dynamical forming of 
a strongly entangled state of bosons. 



II. SEMICLASSICAL VERSUS QUANTUM 
DYNAMICS OF TWO COUPLED INTERACTING 
BOSE SYSTEMS 



junction described by a free harmonic oscillator. The in- 
teraction A is responsible for the anharmonicity. Note 
that for A < 1 the solutions n — 0, = 0, 7r are sta- 
tionary; i.e. the phase difference between the two wells 
can be either or ir. On the other hand— for A > 1 the 
solution with — n becomes unstableEStj, and instead 
the new minima appear at 
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We will now consider the properties of the two well 
system for the two classes of initial conditions discussed 
in Section I in turn. Each subsection below also contains 
a comparison with the exact results obtained by a full 
quantum solution of H. 



A. Mott insulating initial state 



The comparison between the semiclassical and quan- 
tum theory of the two-jaiell system has been presented 
earlier by Milburn et alEB, although for initial conditions 
different from those we shall consider here. 

First we will focus on the semiclassical description of 
the two well system, when the total number of bosons is 
much greater than 1. In this case the Gross-Pitaevskii 



equations implied by (1.4) and (1.5) are 



l ir = -^2 + A|^ii 2 V'i 

at 



(2-1) 
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The total number of bosons |"0i| 2 + I 1 2 is a constant 
of the mot ion; with our normalization for ipj described 
above (|T^) , we have l^il 2 + \^\ 2 = 2. 
We use the parameterization: 



^1,2 = VlT^e ie ^ 



(2.3) 



Note that only the relati ve p hase o f yb\ and t />2 i s an 
observable. Substituting (2J3) into (2T) and ( |2.2| ) we 
obtain: 
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+ An + AXn \J\ — n 2 cos = 0, 
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After further manipulation this system reduces to a 
single second order differential equation for the continu- 
ous variable n: 



<Pn 
dt 2 



+ An + 4An cos 



Xn 2 







(2.6) 



with initial conditions: n(0) = rip^ini^S) / 'dt — 2sin0 o . 
Similar equations were derived int2lEl Without interac- 
tion (A = 0) we have a situation of a single Josephson 



As in Section I A , let us assume that initially the two 



condensates are completely uncoupled. We will consider 
their evolution in the semiclassical and quantum calcula- 
tions in turn: 



1. Semiclassical theory 



;From the discussion in Section I A , we have Uq = 



and 0o is a uniform random variable. We will study the 
correlation between %p\ and i\>i as a function of time. It 
is easy to show that 



(2.8) 



where the average is taken over all possible initial phases 
0o- The correlator is proportional to the product of the 
coupling constant A and the variance of n, reflecting the 
usual phase-number uncertainty relation. 

Before proceeding with quantitative analysis let us ar- 
gue qualitatively what happens with the system. Sup- 
pose A< 1. Then ( |2.6D is equivalent to the motion of a 
particle in a harmonic potential with random initial ve- 
locity. Because the frequency of the harmonic oscillator 
doesn't depend on the amplitude, (n 2 (t)} is a periodic 
function of time wit h T = ir/2. If A is still small but 



not negligible, then (2.E) still describes motion in a har- 
monic potential, which, however, depends on the initial 
conditions. As a result the oscillations of (n 2 (t)) become 
quasiperiodic and decay with time. In the limit of large 
A the oscillations completely disappear and the steady 
state solution develops during the time t ~ l/\/A- 

For weak coupling A, equation ( |2.6D can be solved ex- 
plicitly. Thus for A = 



(n 2 (t)) 



1 - cos At 



(2.9) 



For small A the approximate analytical solution is: 
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1 + A cos < 



(2.10) 



It is easy to see that at large t we have the following 
asymptotic behavior: 
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(2.11) 



so that the variance of n approaches the steady state 
value of one fourth. We note that the amplitude of oscil- 
lations decays with time as t~ 1 / 2 and on top of that there 
are beats with the characteristic frequency LObeats « 4A 
(see Fig. |]). For large A the oscillations decay very 
rapidly and (n 2 (t)} quickly saturates at the steady state 
value, which decreases with A (see Fig.ll|). 



2. Quantum theory 

Let us now study the quantum case. The Heisenberg 
equations of motion are: 



(2.12) 



where square brackets denote c omm utator, j = 1, 2 and 
the Hamiltonian H. is given by (1.1). It turns to be con- 
venient to use the following Heisenberg operators: 



$ = a\ai — a\a2, 



2'ia + &\ci2j 

a>2<i2 — a\ai. 



(2.13) 



We introduce hats over the operators to distinguish them 
from numbers appearing in the semiclassical treatment 
and expectation values of the operators . It is easy to see 
that the following combination 



* n l = * rr 

2N 4J 



(2.14) 



commutes with the Hamiltonian. Using this fact the sys- 
tem fl2.12j ) can be reduced to a single differential equa- 
tion: 



cPn 
~a¥~ 



Ah ■ 



2A 



N 2 



(2n 3 - { 



n,h 2 s) + 



with the initial conditions 
n(0) = h s , 



dh 
HI 



= -2i$ s 



} + ) = o 

(2.15) 



(2.16) 




A 0.15 



0.10 



0.05 



0.00 



40 60 

Time 



100 




10 12 




FIG. 1: Semiclassical variance of n as a function of time. The 
insert on the top graph has a different time scale. 



In the equations above {...}. denotes the anticom- 
mutator, and the subindex s means time-independent 
Schrodi nger operators. We note that the second rela- 
tion in ( 2.16 ) holds for all times if we use $ instead of 



t=o 



In the noninteracting case (A = 0) the solution of (2.15) 
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h(t) = n s cos 2t — i$ s sin 2t. 



(2.17) 



The initial conditions corresponding to the ground state 
for A > X S i is \I) = \N/2,N/2). Note that such a state 
is possible only if N is even. The generalization for N 
odd is straightforward, but we will not do it here, since 
our major goal is to compare quantum and semiclassical 
pictures. Simple computation shows that 



n 2 (t) 
N 2 



1 

7V2 



(I\n 2 (t)\I) = 



l-cos4iiV + 2 
4 N ' 



(2.18) 



Comparing ( 2.18 ) and (2.9) we see that the only differ- 
ence between the semiclassical and quantum results in 
the noninteracting case is the presence of an extra nu- 
merical factor 1 + 2/N in fl2.18|). 



J 



In the weakly interacting regime (A < C 1) we can ne- 
glect terms proportional to A 2 . Then ( [2.15 ) simplifies 
to: 



d 2 h A „ 2A 



{n,* s } = 0. 



(2.19) 



It is very convenient to solve this equation in the eigen- 
basis of 



\k) 



2-JV/2 



y/k\(N-k)\ 



(«L + 4) fc («L-4)^ fe |o), 



(2.20) 



where k = 0,1,. ..N. One can show that for the initial 
Fock state \I) — \N/2, N/2) the variance of n is: 
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(2.21) 



Comparing ( [2.21 ) and ( [2.11 ) we see that in contrast to 
the continuous integral in the semiclassical case there is 
a discrete sum in the quantum. One can formally obtain 
( p. 11 ) from (2.21) in the limit N — > oo using Stirling's 
formula and transforming the summation over k to inte- 
gration. It turns out to be more convenient to normalize 
the variance of n to N(N + 2) instead of N 2 . If the to- 
tal num ber of particles N = 2, there is only one term 
in ( 2.21 ), so the oscillations are completely undamped. 
For N = 4, there are two terms and we expect perfect 
beats; i.e. the amplitude of oscillations first goes to zero 
then completely restores and so on. For N > 6 there 
are several terms contributing to the sum. At relatively 
small time scale X 2 t/N <C 1 frequencies in different terms 
are approximately equidistant: Ail w 8X/N so the am- 
plitude of oscillations is a periodic function. However at 
a larger time scale the phases become random and peri- 
odicity disappears. Figure |2|(a) shows the comparison of 
the variance of n for N — 2 and N — 4 with the semi- 
classical result. On short time scales already N — 4 gives 
an excellent agreement. In fact the semiclassical and the 
quantum curve (for N = 4) are completely indistinguish- 
able. The behavior of the amplitude of oscillations of 
n 2 is plotted in Fig. ||(b). It is clear that with increas- 
ing TV, the semiclassical approximation works for longer 
and longer time scales (see also Ref. |23|). However in 
a quantum system the recurrence time is always finite, 
so ultimately at t > 1/Afi, the semiclassical description 
breaks down. 

In Fig. H we present the numerical solution for the case 
of intermediate and strong couplings. As was discussed 



before for small N, the amplitude of oscillations fluctu- 
ates, being completely chaotic at large time scales. How- 
ever, at sufficiently small time, the oscillations gradually 
decay, approaching the semiclassical result. At interme- 
diate times the amplitude of the oscillations experiences 
beats (compare with Fig. ||). Note that for the large 
coupling, the semiclassical description breaks down very 
early. 



B. Modulated phase initial state 



tion 



W e tu rn next to the initial conditions described in Sec- 
[B, where the initial state has a phase order. In 
semiclassical picture n and 4> are commuting variables 
and we can fix them at t = independe ntly . For simplic- 



ity let us consider uq — 0. Then from (2.6) it is obvious 
that only (f> — 0, it give the stationary solutions. As 
we discussed above, n — and <j> = is automatically 
a ground state for all positive values of interaction A, 
therefore it is always stable under small fluctuations. On 
the other hand if 4>q — it then n = is (meta)stable for 
A < 1 and unstable for A > 1 (see Ref. [l(] for the de- 
tails). Suppose that we start from = 7r, n = Q, A = 
and adiabatically increase A. Then n 2 remains close to 
zero while A remains smaller than critical value. After 
that n 2 rapidly increases and the system spontaneously 
goes to the Schrodinger cat state, where all the bosons 
are either in the left or in the right well. A similar pic- 
ture holds in the quantum mechanical description. The 
principal difference is that instead of a sharp transition 
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FIG. 2: (a) Semiclassical (solid line) and quantum variance of 
n as a function of time for the weak coupling case A = 0.05. 
Dash line corresponds to the total number of particles N = 2, 
dot line does to N = 4. Solid and dot line are indistinguish- 
able on this plot, (b) Amplitude of the oscillations of the 
variance of n versus time. 



FIG. 3: Variance of n as a function of time for intermediate 
(a) and large (b) coupling constants. Note that for larger N 
semiclassical approximation works well for longer time scale, 
but eventually always breaks down. 



at A = A c , there is a smooth crossover between the ini- 
tial and the final states. Fig. |i| shows the variance of n 
as a function of time. For comparison we consider both 
symmetric (0 = 0) and antisymmetric (c/> = n) initial 
conditions. 



The equilibrium number of bosons in the central well 
(j = 0) is N, and so \ipo\ 2 = 1 in the Mott insulating 
ground state. 

We divide our discussion according to the initial con- 
ditions considered in Section I. 



III. SEMICLASSICAL DESCRIPTION OF 
MULTI-WELL BOSE GASES 



The full quantum solution of the many well case 
rapidly becomes numerically prohibitive with increasing 
N, and so we will confine our discussion in t his se ction 
to the semiclassical GP equation. From ( |l.4| ) and (1.5) 
this is 



■A 
' at 



(^•+1 + ty-l) + ^jrTpj + A|Vi|Vi.(3.1) 



A. Mott insulating initial state 



We will compute the correlation function V g (t) defined 
in ( |l.9| ) for two limiting possibilities for the weight func- 
tion g: g(j) = dj^i and g(j) = const, where in the former 
(latter) case one computes the nearest neighbor (global) 
phase correlation. Using the GP equations Q3 . If) we can 



8 



A 0.6 

N 

c 

V 



a ) 


A,-U.UUO t — 


N=2 

N=4 

N 16 

N=64 ,'/ 

GP "7 

i'i : 
•i I 




:i : 

— ^ -*fi/ 

„--' •'/'' 


Antisymmetric 
Initial State 



0.25 



A 0.15 



C 

V 



0.05 



0.00 



200 400 600 800 1000 

Time 



Symmetric Initial State 
^=0.005 t 



•N=2 
■N=4 
N=16 
• N=64 



200 400 600 800 1000 

Time 



1. No external potential and periodic boundary conditions 

Let us assume that the lattice forms a periodic ar- 
ray of quantum wells and there is no external potential 
(Vj = 0). For the nearest neighbor correlation similarly 
to the two well case it is easy to show that 



A 



i) 2 . 

(3.3) 



This equation shows that the nearest neighbor coher- 
ence is proportional to the product of the coupling con- 
stant and sum of the variances of number of bosons in 
each well. From the previous section we can expect that 
if the interaction is weak, then variances of rij at short 
time scales will be fluctuating and governed by the nonin- 
teracting tunnelling Hamiltonian. With increasing time 
the interaction will suppress the fluctuations leading to 
some steady state. In the noninteracting case, (3.1) is 



just an ordinary Schrodinger equation, with eigenstates 



a 27Tikj/N 



M 



corresponding to the eigenenergies 

27rfc 



E k 



-2 cos ■ 



M 



(3.4) 



(3.5) 



FIG. 4: Variance of n for the two wells for adiabatically in- 
creasing interaction \(t). The initial state is (a) antisymmet- 
ric (4> = 7r) and (b) symmetric (<fi = 0). 



Here M is the number of wells. Expanding the initial 
insulating state in terms of the eigenstates defined above 
and propagating them in time we obtain 



show that 

dV g {t) 
dt 



J ^(V,- + A|^W| 2 ) ff (|j-£|) 



(3.2) 



Note that for uniform potential T> g (t) changes only due 
to the interaction. In this case, the ratio T> g (t)/\ has 
a hnite limit at A — ► 0. We will consider the solution 
for T> g (t) with and without an external potential in the 
following subsections. 



jV 



J2(\Mt)\ 2 - l) 2 = M \ 1 - ]T \F(j, t)| 4 ] , (3.6) 



where 



N-l 



kj/M+t cos 2-Kk/M) 



k=0 



(3.7) 



For several different values of M the function Dg(t) 
at vanishing A is: 
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vf(t) = 


A . 2n 

- sin 2 2t, 
2 




(3.8) 


vl{t) = 


OA ,„ „ . . 9 O 

— (2 + cos3t) sin 2 -r, 
9 2 




(3.9) 


v%{t) = 


^(7 + cos2i) sin 2 2t, 




(3.10) 




4A nn o /?, /=■. „ 5^ 3\/5 
— (10 — 2 cos v or — cos v 5r — 2 cos -t cos — — 
25 2 2 


cos V 5c), 


(3.11) 


v & g {t) = 


A , 

— (63 - 8 cos t - 12 cos 2t - 24 cos 3i - 6 cos At - 
36 


- 12cos6t — cos8t), 


(3.12) 


V g M (t)- 


"^( i_jo(<)4_2 £ Jm( * )4 ) at m ^ 


oo. 


(3.13) 



Clearly T>g{t) is a periodic function only for M — 
2,3,4,6 (this is, in fact true, not only for the nearest 
neighbor case). For many wells the number of harmonics 
contributing to the variance of n becomes large and os- 
cillations become more chaotic and weaker in amplitude. 
In the limit M — > oo, T) g l (t) is a monotonically increas- 
ing function. If we add the interaction, then the overall 



picture remains similar to the two well case. Namely, 
for small A the amplitude of oscillations slowly decays in 
time. For strong interaction, the variance of n reaches 
steady state value in a very short time scale. 

In the opposite to nearest neighbors limit g(\j — £\) = 
const, one can show that at A — > 



For example 



M 2A y-T sin 2 r(l + cos(27rfc/M)-cos(27rm/Af)-cos(27r(fc-m)/Af)) 

9 y> ~> m 1 + cos(27rfc/M) - cos(27rm/M) - cos(27r(fc - m)/M) ' ^ ' ' 

K> Tit — 



£>*(*) = ^ sin 2 2t, (3.15) 



■" ' 2 
A 

45' 
X_ 

160' 
1 

240' 



D|(t) = _ (3-2cos3t-cos6r), (3.16) 

V*(t) = -^(13- 12cos4i-cos8£), (3.17) 

T^ait) = A-^- (33 + 16 cos i - 24 cos 2t - 8 cos 6t - cos 8i) , (3.18) 
w 240 

V?® - ™ r f\e ld e^ 2 ^ + CO ; e ^ C °; 9 ^ C ^-J^ at M ^ oo. (3.19) 



The behavior of T> g (t) at large M is very different 
for nearest neighbor and global correlations (see Fig. ||). 
While the former rapidly reaches a steady state value, 
t he la tter oscillates in time. Indeed the denominator in 
( |3 .14 ) selects only low frequency harmonics in T> g , freez- 
ing out high frequency oscillations, especially at longer 
time scales. 

Figs. ^| and [j] show T> g {t) for six and twelve wells re- 
spectively. Six wells give periodic time dependence, while 
N = 12 corresponds to chaotic behavior. Note that in all 



cases high frequency modes are suppressed for the case 
of global phase correlations. 



2. Parabolic confining potential 

So far, we have considered the rather hypothetical sit- 
uation of quantum wells sitting on a ring. However, 
usually one achieves confinement using a trap, which is 
equivalent to a nonuniform external potential Vj in (|3.1|). 
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FIG. 5: Time dependence of the coherence T> g (t) for the 
weakly interacting Bose gases at large number of wells (M — > 
oo). Note that nearest neighbor correlation rapidly saturates, 
while the global coherence exhibits oscillations. 



The most common shape of this potential is parabolic 
(Vj oc j 2 ) and we focus on this case, although the anal- 
ysis of other potentials is similar and straightforward. 
As before, we will first study the non-interacting system: 
(A = 0). 



.dipj 



if 



(3.20) 



This is a linear Schrodinger equation with stationary 
states found from 



Eipj = -(tp j+ i + i/>j. 



(3.21) 



In the Fourier space the same equation looks more famil- 
iar: 



Eip(k) = -2cosfc^(fc) - 



2 dk 2 ' 



(3.22) 



describing the motion of an one-dimensional particle of 
mass living on a circle with the external potential 
U(k) = — 2cos(fc). Note that the same type of equation 
describes Josephson junctions with charging energy. If 
the parabolicity is weak (£ <C 1), then the bosons form 
closely spaced extended states at low energies. In the 
Fourier space this is equivalent to having a heavy parti- 
cle in the —2 cos A; potential. With a good accuracy one 
can describe the energy spectrum inside such a well using 
the WKB approximation. This is justified both for low 
energies, where — 2cosfc m — 2 + k 2 and the WKB gives 
the exact energy spectrum and for high energies WKB 
works well for any potential. In fact there is a little sub- 
tlety near energy close to 2, since the potential there is 
almost flat and can not be approximated by a linear func- 
tion, but this is not very important. So the approximate 
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FIG. 6: Time dependence of T> g (t) for 6 wells; solid and dash 
lines correspond to nearest neighbor and global correlations 
respectively. Without interaction (A — > 0) T> g (t) shows reg- 
ular periodic behavior in time. Nonzero interactions leads 
to decay of oscillations. High frequency oscillations of global 
correlation function are effectively suppressed. 
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FIG. 7: Same as in Figure (g) but for 12 wells. Without 
interaction oscillations are chaotic. Low frequency dominate 
the global correlation function here as well. 
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WKB spectrum is given by 

-cos -1 E/2 



— 7T + COS 



E/2 



(E + 2cos/c) dk = ir(n + 1/2) 



(E + 2 cos k) dk = 2nn, 



(3.23) 



where the top (bottom) equation corresponds to E < 2 
(E > 2). In the first equation even or odd n de- 
scribes even and odd states (in both real and recipro- 
cal space), respectively. For energies E > 2, the second 
equation gives complete degeneracy between even and 
odd energy levels. In real space roughly all states with 
E > 2 are localized in individual wells, and degener- 
ate while those with E < 2 are spread through many 
wells. Fig. ||(a) briefly summarizes this discussion show- 
ing the exact spectrum for £ = 0.1 (The WKB result is 
indistinguishable by eye from this graph). Clearly the 
low energy levels are approximately equally spaced, re- 
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vealing the famous property of a harmonic potential, the 
spacing decreases as the energy approaches 2, and starts 
linearly increasing for E > 2 as in a usual square well. 
If £ > 1 then bosons become localized within individual 
wells and their energies follow external potential. The 
crossover from weak to strong parabolicity is a finite sys- 
tem analog of the Anderson transition. It is important 
to note that this is a purely semiclassical transition in 
this case, because it is derived in the Gross-Pitaevskii 
picture. The "quantum mechanics" here originates from 
the wave nature of the classical field ip. If the average 
number of bosons per well is much larger than one, then 
the semiclassical picture, where number of bosons and 



their phase commute, holds until the typical fluctuations 
of ip 2 becomes of the order of 1/N -C 1. This occurs 
deep inside the insulating regime, where the energy in 
GP approach is anyway almost phase independent. 

After deriving the energy spectrum we can proceed 
with study of the dynamics of the condensate. Note that 
( |3.2| ) yields that time derivative of T> g (t) is not equal 
to zero even without interaction (A = 0). Therefore we 
anticipate that the results for the parabolic and flat po- 
tentials will be strongly different, at least in the weakly 
interacting regime. If the initial phases are uncorrelated 
then it is not hard to show that at A = 



V g (t) = 2j2v(j)g(\3 - £\) £ N^l(j)^a(p)^(p)M^ 



. 2 Ea-E a , 

sin — ^ — t 



E a 



(3.24) 



where Nq is the initial number of Bosons in the well 
number p, ip a and E a are the eigenfunction and energy 
of the level a respectively . If starting from the ground 
insulating state then 



N$ = 1 - for V p < n, 
A* 

A P = for V p > n, 



(3.25) 
(3.26) 



with [i being a che mical potential. Let us make few com- 
ments about ( 3.24 ). Levels (3 and a must have the same 
parity, meaning the lowest harmonic contributing to the 



sum will be w„ 



2min a (E a+ 2 — E a ) > 0. Because 



Nq is centered near the bottom of the well, only lev- 
els with delocalized wavefunctions will contribute to the 
sum. In particular, degenerate levels with E > 2 can 
be safely thrown away. If g(\j — £\) is constant, then 
summation over m ensures that the major contribution 
comes from = 0; therefore T> g (t) contains mostly har- 
monics with lu — E2 ~ Eq, lu — E4 — Eq, etc., with the 
strongest weight at the smallest frequency. Note that 
at small energies and weak parabolicity the lowest en- 
ergy levels are approximately equally spaced, therefore 
the whole expression for T> g (t) will be a quasi-periodic 
function of a frequency lu ss E2 — Eq. However, because 
this equidistance is not exact, the periodicity will be only 
approximate, and at a short time scale the amplitude of 
oscillations will slowly decay. On the contrary for the 
nearest neighbor phase coherence g(\j — £\) = 8j t t±i nei- 
ther (3 nor a are bounded to the ground state and we 
expect that all kinds of allowed frequencies E a — Ep will 
give contributions. Clearly in this case dephasing occurs 
much earlier and the amplitude of oscillations is much 
weaker. Also the characteristic frequency of the oscil- 
lations for the nearest neighbor case will be somewhat 
larger than that for the global case since the level sepa- 



ration decreases with energy. Fig. ^ shows time depen- 
dence of V g for nearest neighbor and global correlations 
at the parabolicity £ = 0.08. From the above analysis we 
should expect the major oscillations at the period 



2?r 



Eo. — En 



2£ 



8, 



(3.27) 



which is indeed very close to the numerical value. 

Interesting things happen if we turn on the interaction. 
In particular, if A is of the order of one, the oscillations 
become much more pronounced and smooth compared 
to noninteracting case (see Fig. ||). This is at first quite 
an unexpected result, since we know that the interaction 
leads to decoherence and saturation ofT> g . However this 
is not the whole story. In the previous analysis we saw 
that at least for the T> g (t) : interaction "kills" high fre- 
quency contributions first. But that is precisely what we 
need for harmonic behavior. So crudely speaking, small 
or intermediate interaction removes harmonics causing 
dephasing of the noninteracting function T> g . If interac- 
tions become strong A > 1, then the noninteracting pic- 
ture is irrelevant and we come back to the usual behavior 
with fast saturation of V g . Notice from Fig |[ that the 
noninteracting and interacting pictures are quite different 
at small time. This can be also understood naturally as a 
result of interplay of many harmonics at early stage of the 
evolution. Hence we expect that the typical time scale 
for the first maximum in the interacting problem will be 
of the order of the tunnelling time, which is much shorter 
than inverse level spacing. However at later times only 
slow harmonics survive leading to slight modifications of 
the noninteracting picture. 
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FIG. 9: Time dependence of T> g for the nearest neighbor cor- 
relation (a) and global correlation (b). The period of oscilla- 
tions scales as l/\/£ an d the amplitude is finite even without 
interaction A = 0. At large A T> g (t) saturates very fast sim- 
ilarly to the flat potential. At intermediate coupling A ~ 1, 
however, the oscillations become more pronounced than in 
the noninteracting regime. 



B. Modulated phase initial state 

It is also s traightforward to generalize the discussion of 
Section II B to the case of the periodic lattice. Namely, if 
the number of wells is even, then the state with a relative 
phase shift 7r, and equal numbers of bosons in the wells, is 
metastable for weak interaction. If A increases gradually, 
then when it reaches a critical value A c , this state be- 
comes unstableoEi The critical value of A can he,fAund 



from the linear analysis of (3.1) near the 7r stafr 



v e 



iqj — icjt\ 



(3.28) 



where u and v are the small amplitudes and q =/= is 
the wave vector of the perturbation. Substitution of this 



for the eigenfrequencies u>: 



lo + 2 — 2 cos q — A 
-A 



which has two solutions 



-A 

2—2 cos q — A 



= ±2V(1 - cosg) 2 - A(l - cosg) 



0, 

(3.29) 
(3.30) 



Clearly w is real if A < 1 — cos q. Otherwise, fluctuations 
with wavevector q become unstable since the frequency 
becomes complex. The lowest nonzero q for the periodic 
boundary conditions is 2tt/M, so the critical value of the 
interaction, where the it state becomes the saddle point 
rather than local minimum is 



X r = 2 sin 



M 



(3.31) 



expansion into (3.1) gives the following secular equation 



Similar to the two well case, the bosons undergo a spon- 
taneous transition to the superposition of states, where 
all of them are in one of the wells. The time dependence 
of the variance of N is analogous to that plotted on the 
top graph of Fig [| (see Fig |l0|). We remark that a "slow" 
or adiabatic increase of interaction must be understood 
carefully. In the GP picture, an adiabatic increase of in- 
teraction means that the characteristic time scale is much 
smaller than the tunnelling time: (din X/dt <C 1). On the 
other hand, for the quantum problem adiabaticity would 
imply that din X/dt is much smaller than the level spac- 
ing, which is proportional to inverse number of bosons. If 
the interaction is increased adiabatically in the quantum 
mechanical sense, then the system would follow the local 
minimum of the metastable state, and when A becomes 
larger than the critical value, it will undergo a sponta- 
neous transition to the dipole state (or a superposition of 
the dipole states) with broken translational symmetry. 



IV. CONCLUSIONS 

We have studied the non-equilibrium temporal behav- 
ior of coupled bosons in a lattice. We predicted dynam- 
ical restoration of the phase coherence after a sudden 
increase of the tunnelling in a system initially in a Mott 
insulating state. In the strongly interacting case, A> 1, 
the coherence reaches a steady state rapidly (within a 
Josephson time). On the other hand, time evolution in 
the weakly interacting regime A < 1 depends strongly on 
the details of the confining potential. We predicted that 
in a parabolic potential Vj = £j 2 /2 the coherence exer ts 
decaying oscillations with period T cx l/\/f (see ( 3.27 )). 
The period and the amplitude of oscillations only depend 
weakly on interaction in this case. On the other hand, 
if the confining potential is flat, then the oscillations are 
either periodic (for a particular number of wells in a lat- 
tice) or chaotic. Here the interaction leads to the decay 
of the oscillations with time. In both cases the system 
ultimately reaches steady state with nonzero coherence 
(dynamical Bose Einstein condensate). 
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FIG. 10: Sum of the squares of number of bosons in different 
lattice cites (with normalization ^ . rij = 1). Clearly uniform 
distributi on is stable until interaction is smaller than the crit- 
oo we have £V i "j ~~ > 1 implying 



ical value 3.31. At t 



that all the bosons populate one of the wells. 



For the two well case we explicitly tested the validity of 
GP approach. It was shown that the mapping of the de- 
terministic quantum mechanical motion to the stochas- 
tic GP equations is essentially exact for time less than 
the characteristic inverse level spacing t < N/X. Apart 
from the slight renormalization of the overall constant, 
the mapping is excellent in this time domain already for 
two bosons per well. For stronger interactions, the semi- 
classical and quantum mechanical trajectories start to 
depart faster, as expected. 

We also considered the dynamical appearance of 
"Schrodinger cat" state under a slow increase of inter- 
action from an initial phase modulated tt state. The tt 
state is stable while interaction is weak and becomes un- 
stable when A > A c . In the GP picture, this instability 
leads to the symmetry breaking, so that all the bosons 
spontaneously populate one of the wells. Quantum me- 
chanically this means that the final configuration is the 
superposition of states in which bosons occupy different 
lattice sites. This approach can be used experimentally 
for the creation of strongly entangled states. 
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APPENDIX A: MEAN FIELD GROUND STATE 
OF THE BOSON LATTICE SYSTEM IN A 
PARABOLIC POTENTIAL 



The problem of the Mott insulator transitions for infi- 
nite arrays of bosons have been extensively studied duiu 
ing the last decade, see for exampleBotl It was showncl 
that the mean filed calculations qualitatively captures 
the two possible phases and gives a good estimate for the 
phase boundary. Recently, using quantum Monte-Carlo 
methods, an exact ground state for— the system of bosons 
in a parabolic potential was founoH It was shown that 
near the expected transition, the global compressibility 
does not vanish due to the spatial inhomogeneity. How- 
ever, still the bosons form local insulating domains sepa- 
rated by narrow superfluid regions. The Monte Carlo ap- 
proach, though very powerful, is incapable to solving the 
problem with many bosons per well. Therefore we think 
that for qualitative understanding of the ground state as 
a function of the interaction strength, it is worthwhile to 
do a mean field calculation. 

The details of the derivation of the mean field equa- 
tions can be found in Ref. ^. Here we will only outline 
the principal steps. 

The mean field version of the free energy, correspond- 
ing to (1.1 ) is 



+—ala j (a i j a j - 1), 



(Al) 



where fi is the chemical potential. The variational pa- 
rameter bj, corresponding to the ground state is: 



W+i 



(A2) 



where the average is taken in the ground state of (Al) 
We can define the order parameter 



P = J2 b Pr 



(A3) 



The self consistent evaluation of the mean field bj is 
straightforward and the resulting order parameter is plot- 
ted in Fig [ll]. The graph (a) corresponds to few bosons 
per lattice site. If the interaction (U) is strong enough, 
then the order parameter, forms a domain structure sim- 
ilar to that predicted inE3. For a large number of bosons 
per well, the quantum fluctuations start playing a role 
when U becomes of the order of the number of bosons in 
the central well (N sa fi/U), and the smooth GP shape 
of boson density (p) breaks down. For very strong inter- 
action, the actual profile of p becomes sensitive to small 
variations of the mean density of bosons per central well. 
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FIG. 11: Mean field order parameter for different interactions 
in a parabolic potential (Vj = £j 2 /2). Graph (a) corresponds 
to few bosons per site and the other two graphs do to many 
bosons. 



APPENDIX B: AMPLITUDE FLUCTUATIONS 
NEAR THE SUPERFLUID-INSULATOR 
TRANSITION 

This appendix reviews results on the damping of the 
amplitude oscillation mode near the superfluid-insulator 
transitioru-JTLOtivated by the recent paper o f Al tman and 
AuerbachHl. As we discussed in Section [A, we have 



considered a system deep in the Mott insulating phase 
(with A S> As/) taken suddenly to parameters for which 
the ground state was deep in the superfluid phase (with 
A <C As/), while Altman and Auerbach consider the case 
when both the initial and final values of A were not too 
far from As/, but remained on opposite sides of it. 

A key ingredient in the dynamics of the amplitude 
mode for A < As/ is the damping induced by emission 
of the Goldstone "spin wave" or "phonon" modes. This 
problem was considered in Refs. |l8| , p6| , and it was found 
that the amplitude oscillations were overdamped in the 
A < As/ scaling limit associated with the second-order 
superfluid-insulator transition. We will review these re- 
sults below, and display expressions which also allow us 
to move beyond the scaling limit to values of A much 
smaller than As/ (see [fjg] ); the amplitude mode can 
become oscillatory in the latter regimetlrEa. This also 
consistent with the considerations of the present paper, 
where we have found that the oscillations of the super- 
fluid coherence were present in the parabolic multiwell 
case for A = 5 in Fig ^|, but were fully overdamped for 
A = 10 (not shown). We found similar behavior in the 
complete quantum solution for the two-well problem — 
however in the latter case, the oscillations reappeared at 
very large A ~ N 2 : these are the "number" oscillations of 
the Mott insulator, and were also found in Ref. [l8| The 
fate of these very small and very large A oscillations in the 
multiwell case near As/ requires a treatment of the inter- 
acting quantum dynamics: this was done in Refs. fl8|]2(i| , 
and the results are reviewed here. 

As is well known, we can describe the superfluid- 
insulator transition by the N = 2 case of the N- 
component ip field theory, where the superfluid order 
parameter ip in Section I 



ip ~ (pi + i<f2 
The action for A close to As/ is 



(Bl) 



S = 



d d xdr 



1 (Y7 \2 , 1 fa ^2 ( r c + s) 2 
^{V x ip a ) +^{dr^ a ) g fa 



2N 



(B2) 



where a = 1 . . . N, c is a velocity, d is the spatial dimen- 
sionality and u is a quartic non-linearity. The coefficient 
of ip a is used to tune the system across the transition, 
and the value of r c is chosen to that the transition oc- 
curs at s = i.e. s ~ A — As/- We assume that in 
the superfluid phase (<p a ) = NoS at i. The oscillations of 
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the spin-wave modes are given by the transverse suscep- 
tibility x_i_(fc,CL>), while those of the amplitude mode are 
given by the longitudinal susceptibility x\\(k,u)', nere k 
is a wavevector, w is a frequency, and the susceptibilities 
are defined by 



VI (k, w)r ) - N 2 (2n) d+i S(k)6(Lu) (B3) 



Expressions for X±.\\ were given in Refs. |l8j,|26j using 
both perturbation theory in u and the large N expan- 
sion. Here, we collect them with a common notation, 
and interpret them in the present context. To first order 
in u, the position of the critical point is determined by 



2u(N + 2)c 
N 



ld+1 



P 



1 



(27T) 



-l p 2 



(B4) 



where p = (k, —iu)/c) is the (d+l)-dimensional Euclidean 
momentum. In the limit of large N, but u arbitrar y, th e 
value of r c is given simply by the N — > oo limit of 
To first order in u, we obtain for x_l 



2 8csu f 
Xi. (P)=T > - —J 



d d ^ 



1 



{2n) d + 1 q 2 + 2s \ (p + q) 2 



1 



1 

? 
(B5) 



where q is also a (d + l)-dimensional Euclidean momen- 
tum; at N — oo we have simply xj 1 (p) = P 2 ■ The expres- 
sion (B5) describes the spin- wave oscillations, along with 
their essentially negligible damping from their coupling 
to the amplitude mode (as can be veri fied by taking the 
imaginary part of the loop integral in (BE) after analytic 
continuation to real frequencies). 

The damping in the longitudinal modes is much more 
severe, and we will consider it explicitly. To first, order 
in u, we obtain the expression 



Xi\p) 



P 



2s 



4csu(N - 1) 
N 



U(p) + 6xi 1 ( P ). ( B6 ) 



Here the strong damping term has been i nclu ded in Tl(p) 
whose explicit form is discussed below in ( p39| ), while 5\\\ 
contains additional non-singular terms we can safely ne- 
glect. For completeness, we give the expression for the 
latter 



Sxi'ip) = 



\2uc 

N 
36csu 

N 



d d+i 



<1 



{2ix) d + 1 



l 



i 



d d+i, 



2s 
1 



(B7) 



(2^) d +! {q 2 + 2s){{p + q) 2 + 2s) : 



note that these terms always involve coupling to an am- 
plitude mode fluctuation (with "mass" 2s) and this is the 



rea son their contribution is non-singul ar. We find below 



in (B9) that the U(p) contribution in (B6) involves only 
spin-wave fluctuations and hence it becomes very large 
at lo w frequencies, where the perturbative expansion in 
( p36| ) can no longer be trusted. Fortunately, a resumma- 
tion of these singular corrections is provided by the large 
iV expansion, which yields 



Xii 1 (p) =p 2 



2s 



1 + 2cuIi{pY 



(B8) 



it is satisfying to check that (B8) and (B6) are entirely 
consistent with each other in their overlapping limi ts of 
validity of small u and large N. The expression (B8) was 
given earlieilij in the scaling limit, which corresponds to 
ignoring the 1 in the deno minator because n(p) becomes 
large. The utility of (B8) is that it does not have diver- 
gent behavior at small p. 

We turn, finally, to the expression for Tl(p), which is 



n(p) 



d d +h 



2ir) d + 1 q 2 (p 
F d \p\ d - 3 , 



9) 2 



(B9) 



where Fd is a numerical prefactor which is not difficult to 
obtain explicitly. Notice that H(p) is singular as p — > 
in d < 3, and this is the reason for the strong damping of 
the amplitude mode. After analytic continuation to real 
frequences, we have in d — 2 



IL(k,w) 



1 



d = 2: 



(BIO) 



this has a non-zero imaginary part for lu > ck which leads 
to the damping of the amplitude mode. The expression 
for n(p) is infrared divergent in d = 1, and this is the 
signal that there is no true long-range order; nevertheless, 
its imaginary part remains well defined as d \ 1, and we 
find 



1 



4(( W /c) 



(lu — ck) 



d = l 
(BH) 



which again predicts strong d amping at low frequen- 
cies. The expressions (BS-B11) can be used to describe 
the evolution of the weakly damped amplitude mode at 
lu = \/c 2 k 2 + 2s at large s deep in the superfluid, to the 
overdamped mode with no sharp resonance at this fre- 
quency for small s. 
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